Gap inhomogeneities and the density of states in disordered d-wave superconductors 
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We report on a numerical study of disorder effects in 2D d-wave BCS superconductors. We compare 
exact numerical solutions of the Bogoliubov-deGennes (BdG) equations for the density of states 
p{E) with the standard T-matrix approximation. Local suppression of the order parameter near 
impurity sites, which occurs in self-consistent solutions of the BdG equations, leads to apparent 
power law behavior p(E) ~ \E\ a with non-universal a over an energy scale comparable to the single 
impurity resonance energy fio- We show that the novel effects arise from static spatial correlations 
between the order parameter and the impurity distribution. 
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In spite of strong electronic correlations in the nor- 
mal state, the superconducting state of high T c materi- 
als seems to be accurately described by a conventional 
BCS-like phenomenology. The debate over the fc-space 
structure of the BCS order parameter Ak has now been 
resolved in favour of pairing states with d-wave symme- 
try. Since this symmmctry implies that Ak changes sign 
under rotation by 7r/2, there are necessarily points on 
the 2D Fermi surface at which Ak vanishes. The unique 
low-energy properties of high T c superconductors are de- 
termined by the quasiparticle excitations in the vicinity 
of these nodal points. 

For conventional s-wave superconductors, the density 
of states (DOS) p(E) has a well defined gap and is 
largely unaffected by non-magnetic disorder. In contrast, 
p{E) oc \E\ in clean d-wave superconductors, and can 
be substantially altered by disorder. Much of the cur- 
rent understanding of disorder effects comes from per- 
turbative theories, such as the widely-used self-consistent 
T-matrix approximation (SCTMA). In particular, the 
SCTMA is exact in the limit of a single impurity, and 
has been used in studies of the local DOS near an iso- 
lated scatterer jjj. For sufficiently strong scatterers, an 
isolated impurity introduces a pair of resonances at en- 
ergies ±Oo [§,[§, where flo < A is a function of of the 
impurity potential uq and of the band asymmetry. An- 
alytic expressions for Hq(uq) have been given for a sym- 
metric band and in this special instance the unitary 
limit flo — > coincides with uq — > oo. For a realistic 
(asymmetric) band, the relationship is more complex ||. 

For a finite concentration of impurities nt, the SCTMA 
predicts that the impurity resonances broaden, with tails 
which overlap at the Fermi energy, leading to a finite 
residual DOS p(0) The region over which p(E) f» 

p(0), the "impurity band", and has a width comparable 
to the scattering rate 7 at E = 0. In the Born limit, 
the ±f2o resonances are widely separated in energy, and 
the overlap of their tails is exponentially small. In the 
strong-scattering limit, however, the overlap is substan- 
tial and 7 ~ yTArfj where A^ is the magnitude of the 
d-wave gap, T — ni/nN is the scattering rate in the 



normal state, and No is the 2D normal-state DOS at the 
Fermi level. Several recent experiments J9|-|l2j have stud- 
ied quasiparticles in the impurity band in Zn-doped high 
T c materials. Of particular note are attempts to verify a 
provocative prediction Jl^,|ll| that transport coefficients 
on energy scales lo, T < 7 have a universal value, inde- 
pendent of r ftHQ ■ 

The SCTMA which forms the basis of this world-view 
has several limitations. It is an effective medium the- 
ory, in which one solves for the eigenstates of an isolated 
impurity in the presence of a homogeneous mean-field 
representing all other impurities. This approach ignores 
multiple impurity scattering processes which are respon- 
sible for localization physics in metals and may lead to 
novel effects in 2D d-wave superconductors |il||l^ , |i"6[ . 
Another limitation of most SCTMA calculations is the 
use of a 5-function potential as an impurity model. While 
this simplifies the calculation substantially, numerical re- 
sults hint that the detailed structure of the impurity po- 
tential may be important B. A related issue, which 
will be discussed at length in this Letter, is inhomo- 
geneous order-parameter suppression. It is well-known 
that d-wave superconductivity is destroyed locally near a 
strong scatterer, and in the single- impurity limit Jl9| , p0| , 
the additional scattering was found to renormalise Qq 
but, surprisingly, to leave most other details of the scat- 
tered eigenstates unchanged. Here we show, using exact 
numerical solutions of the Bogoliubov-deGennes (BdG) 
equations, that additional novel physics does arise in the 
many-impurity case. 

The main results of this work are summarised in Fig. [I]. 
For a homogeneous order parameter, p(E) saturates at a 
constant value as E — * (in agreement with SCTMA), 
down to a mesoscopic energy scale ~ ^/pL 2 where level 
repulsion across the Fermi surface induces a gap. This 
small gap may be a precursor to a regime associated 
with strong localization in the thermodynamic limit as 
discussed by Senthil et al. Jll|. These authors predicted 
asymptotic power-law behavior in p(E) over an exponen- 
tially small energy scale E2 ~ 1/^(0)^, where p(0) is the 
residual density of states in the impurity band plateau 
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FIG. 1. Density of states for rii — 0.04 and «o = 10. 
Numerical solutions of the BdG equations are shown 
with (BdG+OP) and without (BdG) self-consistent calcu- 
lation of Aij. Inset: T-matrix calculations of p(E) with 
(SCTMA+OP) and without (SCTMA) off-diagonal scatter- 
ing. Also shown is a model with uncorrelated impurity and 
off-diagonal potentials (uncorr.). 



and £l ~ (vf/j) exp(i>i?/i>A + Va/vf) is the quasiparti- 
cle localization length, (vp is the Fermi velocity and va 
is the gradient of A& along the Fermi surface at the gap 
node). The actual value of the DOS at E = is consistent 
with zero but not determined in our work; theoretically 
this point is still controversial. [jl8| 

Figure [l] shows that when the order parameter is de- 
termined self- consistently from the BdG equations p{E) 
is quite different. At low energies the DOS can be fit to a 
power law p(E) ~ \E\ a with nonuniversal a (Fig. |^). The 
power-law is the result of spatial correlations between the 
order-parameter and the impurity potential, and is there- 
fore fundamentally different from those of Nersesyan et 
al. |l5| ] and Senthil et al p6| , where asymptotic power- 
laws were found, with a — 1/7 and a — 1 respectively. 
Unlike the DOS, the dimensionless conductance and in- 
verse participation ratio (Fig. ^) are not changed signif- 
icantly by self-consistency. Finally, we remark that the 

energy scale for the low-energy regime is ?Jf2o, which is 
orders of magnitude larger than Ei for realistic parame- 
ters. 

We employ a one-band lattice model with nearest 
neighbour hopping amplitude t and a nearest neighbour 
attractive interaction V. Substitutional impurities are 
represented by a change in the on-site atomic energy. 
The Hamiltonian is 

(id) a »> CT 
(id) 

where the angle-brackets indicate that site indices i and j 



are nearest neighbours, Ui is the impurity potential which 
takes the value uo at a fraction m of the sites and is zero 
elsewhere, and Ay = — V(cjic^) is the mean-field order 
parameter, determined self-consistently by diagonalizing 
Eq. (Q). Throughout this work, energies are measured 
in units of t, where t is of order 100 meV for high T c 
materials. In non-self-consistent calculations the OP has 
the familiar fc-space form A& = A e i[cos(k x ) — cos(fty)], 
where A^ = \^ ± {Aa± x — ^a± y ) is independent of i. 
Unless otherwise stated, V — —2.3 and p = 1.2 which 
yields Ad = 0.4 (corresponding to vf/va ~ 5) in the 
absence of disorder. Self-consistent solutions show that 
is suppressed within a few lattice constants of each 
strongly-scattering impurity. Throughout this Letter, 
curves marked BdG and BdG+OP refer to the neglect 
or inclusion of self-consistency in Ay . 

The DOS is p(E) = L~ 2 £ Q 8{E-E a ), where L is the 
linear system size and E a are the discrete eigenenergies 
of Ti.. Our numerical calculations were performed on pe- 
riodically continued systems with L < 45. Typical DOS 
curves were obtained for L = 25 by averaging p{E) over 
~ 50 — 500 impurity configurations and ~ 50 — 100 k- 
vectors in the supercell Brillouin-zone. For system sizes 
L > 35, computational constraints restricted us to real 
periodic and anti-periodic boundary conditions. 

An important motivation for the present study is the 
need for a test of the reliability of SCTMA predictions. 
Thus, we use the same disordered lattice model for the 
SCTMA as for the BdG calculations. We would also 
like to model order-parameter suppression within a self- 
consistent T-matrix approximation (SCTMA+OP), Q 
and follow the ansatz of jl!| that the off-diagonal poten- 
tial is 8 A^ = —Aij[8ifl + Sjfl]. This term appears in 
the off-diagonal block of the effective potential. In both 
cases, the T-matrix is a 2 x 2 matrix in particle-hole space 
which satisfies 

Ty (E) = Ua + U iR°( R - R '> E ) T R'] (£)> ( 2a ) 

R,R' 

where G has the Fourier transform, 

G k (E) = [Et - e k r 3 - A kTl - Efc(-E)] ~\ (2b) 

£fc is the tight-binding dispersion, are the Pauli 
matrices, and H^{E) = riiTkk{E). Equation ( |2a| ) is 
solved in real-space to take advantage of the short 
range of the effective potential U. Finally, p(E) — 
— 7r _1 L _2 Im X)fc Gfc(i?), where the trace is over 
particle-hole indices. 

Except for the mesoscopic gap discussed previously, 
the BdG curve in Fig. [j] is quantitatively similar to the 
SCTMA (inset). In contrast, the BdG+OP curve van- 
ishes smoothly as E — * 0, indicating that qualitatively 
new physics has been introduced by the inclusion of or- 
der parameter suppression. To emphasize that this result 
is unexpected, we point to the popular "Swiss cheese" 
model, in which it is assumed that pair-breaking causes 
a pocket of normal metal of radius £q (the coherence 
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FIG. 2. Dependence of p(E) on n, and uq: (a) BdG+OP 
(solid) and SCTMA (dotted) for uo — 5 (top three curves) 
and Mo = 10 (bottom). Also shown (dot-dash) is a fit 
p{E) = y4|i5| a for m — 0.04. (b) Logarithmic plot showing 
power law behavior for rii = 0.04, tto = 5. (c) Dependence of 
power law on wo for m = 0.02. For comparison, SCTMA of 
p(0) vs. uo is shown. Unitary limit is uq ~ 5. (d) Scaling of 
power law with rii. Note that when rii = 0, a = 1. 



length) to form around each impurity. Within this model, 
p(0) is enhanced relative to the SCTMA by ~ uiNqQ. 

We emphasize that the correct spatial correlations be- 
tween an impurity configuration and its self-consistent 
off-diagonal potentials must be preserved for the 
BdG+OP results to arise. This point is illustrated with 
a simple numerical calculation (inset of Fig. [l]) , in which 
Ay is found self-consistently for random impurity dis- 
tributions which are different from those appearing in 
the diagonal block of Ti. The system therefore has two 
distinct types of impurity, one of which is purely off- 
diagonal, with uncorrelated distributions. It is striking 
that there is no hint of the correct low-energy behaviour 
in this calculation. 

It is instructive to compare the BdG+OP result with 
the SCTMA+OP (inset, Fig. 0), since they include off- 
diagonal scattering from the order parameter at different 
levels of approximation. Furthermore, the SCTMA+OP 
does preserve the correlation between impurity location 
and order-parameter suppression. It has been used suc- 
cesfully to describe shifts in VLq due to off-diagonal scat- 
tering Jn|^0| but, here, fails to reproduce the correct low 
energy DOS in the bulk disordered case. Instead, p{E) 
is quantitatively similar to the SCTMA, which is a di- 
rect result of the relative smallness of the off-diagonal 
potential A^/mo « 0.04. 

In Fig. @ we illustrate the dependence of the low en- 
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FIG. 3. Scaling of the inverse participation ratio for 
BdG+OP (solid curves) calculations. A BdG calculation 
(dashed) is also shown for L = 20. Parameters are V = —4.47 
(Ad = 1.34), m = 0.06, and p = 1.0 with between 14 
(L=45) and 500 (L=20) impurity configurations and real 
boundary-conditions. Inset: Scaling inside and outside the 
impurity band. a< is the average of a(E) for \E\ < 0.8. 



ergy DOS on both m and uo- In Fig. ||(a) a series of 
curves shows how the low-energy regime scales towards 
zero-width as m — > for the near-unitary scattering po- 
tential Mo = 5 pH| . That the size of the regime should 
scale faster than 7 with n, is consistent with our earlier 
assertion that the novel behavior stems from a multiple- 
impurity effect. The details of the dilute impurity limit 
depend on the particular value of uq however: at 2% 
impurity concentration, the low-energy regime is signif- 
icantly larger for uq — 10, which lies farther from uni- 
tarity, than for uq = 5. For fixed rii, we find that the 

low-energy regime is ^,£Iq. 

In Fig. ||(b), a logarithmic plot of p(E) reveals that 
the low energy DOS has an apparent power-law depen- 
dence on E, with non-universal exponent < a < 1 
[Fig. ||(c),(d)]. At low impurity concentrations, a is a 
strong function of both and uq, with a a minimum for 
unitary scatterers. For larger n,, a appears to saturate 
at a value which is independent of uq. We assert that 
this power law is fundamentally different from those re- 
ported elsewhere |l5|,[l6|], since it is only observed when 
off-diagonal scattering is present and, as we will see next, 
is unrelated to strong quasiparticle localization. 

We have studied the scaling of both the inverse par- 
ticipation ratio a(E) and the Thouless number g(E) 
with system size. The inverse participation ratio is de- 
fined in the usual way, a(E) — xn(E) / X2{E) 2 , where 
x m {E) = [p{E)L^Y, n A\un{r)\ m + \v n {r)r]8{E-E n ) 
and u n (r) and v n (r) are the particle and hole eigenfunc- 
tions, and is plotted in Fig. for several system sizes. 
a(E) is a direct measure of the spatial extent of the 
wavefunction; it scales as 1/L 2 for extended states and 
is constant for states with localization length £l < L. 
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The scaling and ener gy d ependence of a(E) is in general 
agreement with ref. [|22| where calculations were made 
with similar parameters — there is a crossover in scaling 
between E < 7 and E > 7 consistent with a shorter local- 
ization length in the impurity band. Unlike ref. p2| , wc 
do not find saturation in a< [the average of a(E < 7)] for 
L < 45. For larger disorder concentrations (not shown) 
where < 45 sites, we find that a(E) saturates at E = 
first indicating that £l (E) is an increasing function of E 
in the impurity band. The most important result for 
this work is that a(E), and therefore £l, is not signif- 
icantly different in the BdG and BdG+OP calculations 
despite a substantial difference in p(E). In the calcula- 
tion which is shown, a{E) is actually decreased slightly 
by self-consistency, corresponding to a slight increase in 
the localisation length. 

We have also studied the Thouless number, defined as 
9(E) = E„K(tt) ~ E n {0)]5(E - E n ), where the argu- 
ment of E n refers to the application of periodic or anti- 
periodic boundary conditions in the ^-direction. As L is 
increased, g(E) is expected to cross over from a constant 
in the diffusive regime to exponential scaling indicative 
of strong localization. For L < 45 we find no significant 
scaling of g(E) with L, consistent with what is found 
for a(E). Most significantly, we find that g(E) is nearly 
identical in the BdG and BdG+OP calculations, even 
for low-energy states where substantial changes in p(E) 
occur. 

This behaviour is reminiscent of what is seen in 
Hartree-Fock studies of interacting electrons in disor- 
dered conductors [Q. There, the Coulomb interaction 
enforces spatial correlations between the disorder and 
charge distributions and leads to the formation of a gap 
in p(E) yet leaves the dimensionless conductance 

(a two-particle property related to the Thouless number) 
unchanged. The power-law DOS we observe here may 
have a similar origin; it is certainly clear that p(E) de- 
pends crucially on the spatial correlations between the 
impurity potential and d-wave order parameter. In the 
current work, however, it is the pairing interaction which 
is relevant and the BdG equations provide a mean-field 
description of the pairing interaction which is analogous 
to the Hartree-Fock description of the Coulomb inter- 
action. We speculate that the short-ranged pairing in- 
teraction produces spatial correlations between distant 
impurities via the overlap of the long range tails of 
the single impurity resonances. 

In this work we have shown that spatial correlations 
between order parameter and impurity distributions in 
d-wave superconductors lead to apparent power-laws in 
p(E) at low energies. These results are potentially rele- 
vant to quasi-2D superconductors like BSCCO-2212. Un- 
fortunately, most disorder studies have been performed 
on the anisotropic 3D YBCO system, where correlation 
effects are expected to be less pronounced. Indeed there 
is considerable evidence, particularly in Zn-substitutcd 
YBCO, that disorder does indeed induce a finite DOS 
at the Fermi level [plfKllPj and somewhat weaker evi- 



dence that it scales with disorder in accordance with the 
SCTMA U . Our work should therefore provide a strong 
motivation to study Zn doping, and other types of pla- 
nar disorder, in the quasi-2D BSCCO-2212 system at low 
temperatures. 
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